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A CONTOUR-INTEGRAL BASED METHOD FOR COUNTING THE 
EIGENVALUES INSIDE A REGION IN THE COMPLEX PLANE 


GUOJIAN YIN* 

Abstract. In many applications, the information about the number of eigenvalues inside a 
given region is required. In this paper, we propose a contour-integral based method for this purpose. 
The new method is motivated by two findings. There exist methods for estimating the number of 
eigenvalues inside a region in the complex plane. But our method is able to compute the number 
of eigenvalues inside the given region exactly. An appealing feature of our method is that it can 
integrate with the recently developed contour-integral based eigensolvers to help them detect whether 
all desired eigenvalues are found. Numerical experiments are reported to show the viability of our 
new method. 
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1. Introduction. Consider the generalized eigenvalue problem 

Ax = XBx, (1-1) 

where A, B G The scalars A € C and the associated nonzero vectors x € C" are 

called the eigenvalues and eigenvectors, respectively. Let R be a disk in the complex 
plane and enclosed by circle F. In this paper, we want to develop an approach for 
computing the number of eigenvalues of (EH inside T> exactly. Due to the Mobius 
transformation, the resulting approach can be adapted to the union of intersections 
of arbitrary half plane and (complemented) disks, and so a rather general region [3]. 

In many applications, it is required to know the number of eigenvalues inside a 
prescribed region in the complex plane mui]. For example, it is a prerequisite of 
the eigensolvers based on divide-and-conquer techniques Hi. To get this number, 
the most straightforward way is to compute all eigenvalues inside the target region 
by some method, such as the rational Krylov subspace method m- However, this 
way is always time-consuming and not effective, because we are only interested in the 
number of eigenvalues rather than the eigenvalues themselves. 

When A and B are Hermitian matrices with B being positive definite, i.e., zB —A 
is a definite matrix pencil H; it is well-known that the eigenvalues of EH are real¬ 
valued and lie on the real line [181 HI] ■ Assume that we want to know the number 
of eigenvalues inside interval [a, b]. The standard method works as follows. Compute 
two LDL decompositions: A — aB = LaDaL’^ and A — bB = L^DbL’^. Then the 
difference between the numbers of negative entries in the diagonal of Da and Di, is 
exactly the number of eigenvalues inside [a,b]. The derivation of this method takes 
advantage of the Sylvester law of inertia, see lilllllH] for more details. Obviously, 
the efficiency of this method depends on the accurate computation of two related 
LDL decompositions [Sj. Computing the LDL decomposition requires floating point 
operations of order 0{n^) [TT]; as a result, it becomes computationally prohibitive 
for large-scale problems. When EH) comes to the non-Hermitian case, to the best 
of our knowledge, there is no specihc method for exactly computing the number of 
eigenvalues of (11.11) inside V. 
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The contour-integral based eigenslovers [isiiiniiiiiiii] are recent efforts for com¬ 
puting the eigenvalues inside a prescribed region. The information about the number 
of eigenvalues inside the region of interest is crucial to the practical implementation 
of these eigensolvers. Due to this, recently some methods based on contour-integral 
were proposed to get estimations of this number. Some of this kind of methods make 
use of the stochastic estimation of the trace of spectral operator associated with the 
eigenvalues inside the given region [8l[9l[2l]. However, they may be unreliable in some 
cases, for instance when the matrices A and B are ill-conditioned. In [53], another 
kind of contour-integral based estimate method was presented under the assump¬ 
tion that matrices A and B are Hermitian and B is positive definite. It should be 
pointed out that all these existing contour-integral based methods always just provide 
approximations for the exact number of eigenvalues inside the given region. 

In this paper, we present a contour-integral based method, which can exactly 
compute the number of eigenvalues of (HH) inside V. The derivation of the proposed 
method requires the eigenvalues of (HH) are semi-simple, namely, there are n indepen¬ 
dent eigenvectors, which is always the case in practical situations. Our new method 
is motivated by two findings. The first finding comes from using the Gauss-Legendre 
quadrature rule to approximately compute the spectral operator constructed by a par¬ 
ticular contour integral. More details will be discussed in Section 2. The second one is 
devoted to avoiding the computation of the Weierstrass canonical form of matrix pen¬ 
cil zB — A when using the first finding to count the eigenvalues inside V. We will detail 
the second finding in Section 4. Since our new method is also based on contour inte¬ 
gral, it keeps the promising features of the usual contour-integral based eigensolvers, 
such as having a good potential to be implemented on a high-performance parallel 
architecture. Moreover, it can integrate with the contour-integral based eigensolvers 
dSlEQlEIlEl] to help them determine whether all desired eigenvalues are found when 
they stop. 

The paper is organized as follows. In Section 2, we present the first finding, which 
is derived from using the Gauss-Legendre quadrature rule to approximately compute 
the spectral projector constructed by an integral contour. Since our method needs 
the help of a technique proposed in |24| . we briefly describe the technique in Section 
3. In Section 4, we first detail the second finding, and then give the resulting method 
for counting the eigenvalues inside B. Numerical experiments are reported in Section 
5 to show the viability of our new method. 

Throughout the paper, the following notation and terminology are used. The real 
part of a complex number a is denoted by 5ft(a). We use y/—l to denote the imaginary 
unit. The subspace spanned by the columns of matrix X is denoted by span{X}. The 
rank of X is denoted by rank(X). The algorithms are presented in Matlab style. 

2. Approximate spectral projector. Our discussion starts with the spectral 
projector associated with the eigenvalues inside T, which is constructed by a contour 
integral defined as 



( 2 . 1 ) 


A matrix pencil zB — A is called regular if det {zB — A) is not identically zero for 
allz € C [ZlIIS]. Below is a generalization of the Jordan canonical form to the regular 
matrix pencil case. 

Theorem 2.1 (Weierstrass canonical form [IH])- Let zB — A be a regular matrix 
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pencil of order n. Then there exist nonsingular matrices S and T £ 


TAS = 


Jd 

0 


0 

^n—d 


and 


TBS = 


Id 

0 


0 

Nn-d 


( 2 . 2 ) 


where Jd is a d x d matrix in Jordan canonical form with its diagonal entries corre¬ 
sponding to the eigenvalues of zB — A, Nn-d is an (n — d) x {n — d) nilpotent matrix 
also in Jordan canonical form, and Id denotes the identity matrix of order d. 

Assume that there are n independent eigenvectors, which implies Jd is a diagonal 
matrix and Nn-d is a zero matrix. Let Jd = diag{Ai, A 2 ,..., A^}, with Ai being the 
(finite) eigenvalues [7]. Here the Xi are not necessarily distinct and can be repeated 
according to their multiplicities. 

For z ^ Xi, the matrix {zid — Jd) is invertible. Hence, according to (12.2L the 
resolvent operator 


{zB-A)-^B = S 


= S 


{zid - Jd)-^ 0 

0 {zNn-d - In-d)~^ 

{zId - Jd)-^ 0 

0 {zNn-d - In-d)~^ 

= SD{z)S-^, 


TB 

Id 0 

0 Nn-d 


where 

D{z) = 

and the diagonal block [zid — Jd)~^ is of the form: 


[zld-Jd)-^ O' 
0 0 


[zid - Jd) ^ = 


1 


z — Xi 


z- X 2 


z - Xd . 


(2.3) 


(2.4) 


(2.5) 


Assume that there are s eigenvalues enclosed by F, without loss of generality, let 
them be {Ai,..., Ag}. Then, according to the residue theorem in complex analysis 
[T], it follows from (I2.3D - (I2.5I) that 


Q = S 


27 rv/^ Jr 


® D[z)dz 


s-^ = S 


Is 0 

0 0 


5-1 =5(:,l:g)(5-l)(l:g 


( 2 . 6 ) 


It is easy to verify that = Q, which implies that Q is a spectral projector onto the 
eigenspace span{S'(:_i:s)} corresponding to {Ai, A 2 ,..., Ag} |24) . 

In view of (1^ . the spectral projector Q can be obtained via computing the 
Weierstrass canonical form (cf. ()2.2I) '). However, it is well-known that the Weierstrass 
canonical form is not suitable for numerical computation mm- According to the 
expression (EH), an alternative way is to compute Q numerically by a quadrature 
scheme. In our method, the quadrature scheme is restricted to the Gauss-Legendre 
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quadrature rule [B]. Let c and p be the center and the radius of circle F, respectively. 
Applying the g-point Gauss-Legendre quadrature rule to (12.ip yields 




\ ~ - ^) ^B = S 

j=i 


1 

2 



c)D{zj) 


s-\ 


(2.7) 


Here Zj = c + , Oj = {l+tj)Tr, and tj is the jth Gaussian node with associated 

weight ujj. We remark that Q is an approximate spectral projector. Let 

I 9 

D = - '^ujj{zj - c)D{zj). (2.8) 

i=i 


Comparing to (12.61) and (12.71) , we see that D is an approximation to 
Let p = c + where r £ [0, oo) and 9 G (—tt, tt]. Define 


0 0 


tpip) = 


Jy z — p 


-dz. 


(2.9) 


According to the residue theorem, we know that V’Cm) = 1 when p is located inside 
r, and = 0 when p is located outside F. If ipifJ-) is computed approximately by 
the g-point Gauss-Legendre quadrature rule, then we have 


„ 1 9 I 

tj}{p) ^ V'(m) = - c)—^ 


i=i 






■ 3 - 9 - 

—pcos(tj7r) — V—lpsin(tj7r) 


i=i 


— {pcositjir) + rcos6) — •v/i^(psin(tj7r) -|-rsin0) 


( 2 . 10 ) 


It was shown that '0 (m) ^ ^ if is real-valued and located inside F [23]. This 
observation sheds light on the following theorem. ^ 

Theorem 2.2. If p is enclosed by F, then the real part of ipip) satisfies 

If p is located outside F, we have 


Proof. According to (I2.10p . we have 


Let 


^’{ 9 )] 


1 ^ p^ -I- pr cositj-K — 9) 

2 ^ p^ Pr"^ + 2pr cos(tj7r — 9) 


p^ -I- pr cositj-K — 9) 

p2 _j_ ^2 _|_ 2pr cos(tj7r — 9)' 
4 


( 2 . 11 ) 
















For any given j, one can show that 


^ 1 /O 1 o^ 

2 2 [(p + r cos(tj7r — 0))^ + (r sin(tj7r — 0))2] ^ ^ 

Note that the denominator of the right hand in (12.121) is always positive. It is readily 
to see that 


9jir,e) > 


1 


for r € [0, p) and 9 € (—tt, tt], in which case p is enclosed by F, and 

9j{r,e) < i 


for r € (p, +oo] and 9 € (— tt, tt], in which case p is located outside F. 

On the other hand, it is well-known that = 2 [6l [23]. Therefore, 

^^{ 9 )] = ^J2^i9j{r,9) > i 

i=i 


when p is enclosed by F, and 

^bPi^)] < ^ 

when p is located outside F. □ ^ 

We use Figure Hr] to depict the function 5i['0(p)]. The figure well demonstrates 
the conclusion in Theorem 12.21 

Observe (12.4p . (|2.5p . and (I2.8p - (I2.10|1 . we see that 

-^(2,2) ^ 1 , . . . d. 

Due to this, in the following we always refer to D(j j) as the diagonal entry of D that 
corresponds to eigenvalue Xi,i = 1,... ,d. In viewing of Theorem 12.21 the diagonal 
entries of D can be divided into two categories, that is, > ^}f^i on one 

group and on the other. This fact is the first finding of our work, 

which implies that the number s can be obtained by counting the diagonal entries 
of D whose real parts are larger than i. According to (12.8L getting D requires to 
compute the Weierstrass canonical form of zB — A. However, as was suggested in |2], 
the Weierstrass canonical form is not suitable for numerical computation. 

In our work, we present an alternative method to obtain which does 

not need to compute the Weierstrass canonical form of zB — A. As a result, by 
exploiting the first finding, we can exactly compute the number of eigenvalues of CH) 
inside F. The alternative method is our second finding, since it needs the help of a 
technique proposed in |24j . we briefly introduce the technique in next section before 
presenting the second hnding. 

3. Finding an upper bound for the number of eigenvalues inside F. In 
[24] , a method based on contour integral was proposed for finding a good upper bound 
for the number of eigenvalues inside a given region. Meanwhile, it can produce an 
approximate projection onto the eigenspace corresponding to the eigenvalues inside 
the given region. 
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g = 16 g = 16 



Fig. 2.1. This figure illustrates the function (cf. h2.11 ^ ) when r £ [0,4]. The circle 

r has radius p = 1 with center at the origin. We set the degree of Gauss-Legendre quadrature 
rule q = 16. The left picture shows the general shape of the function, and the right one shows the 
logarithmic scale. 


The method first uses a stochastic estimation of the trace of spectral projector 
Q (cf. ()2.ip i to obtain an estimation of s. By Yp ^ ^nxp{^^ l)i we mean that Yp is 
an n X p random matrix with independent and identically distributed entries drawing 
from standard normal distribution N(0,1). By (|2.6p . one can easily verify that 

^¥.[iT:a.ce{Y* QYp)] = trace(Q) = trace(S'(:_i:s)(5'“^)(i:s_:)) 

= trace((S'“^)(i:s,:)S'(: 4 :s)) = trace(/s) = s. 

Therefore, 

So := \Yace(Y*QYp)'] (3.1) 

provides an estimation for s, see 0111 HU for more details. With this knowledge on 
hand, a method was then given in [24] to seek a good upper bound si of s. Wanting 
to derive the method, we need the following lemma. 

Lemma 3.1 ([24]). Let Y £ If the entries of Y are eontinuous ran¬ 

dom numbers that are independent and identieally distributed (i.i.d.), then the matrix 
{S~^)(^i.m,-.)Y is almost surely nonsingular. 

Let s’l’ be a positive integer and Y),, N„xst(0,1). Consider 

Usi = QYs, = ^-^£{zB-A)-^BdzY,, = 5(:,l:.)(5-^)(l:.,:)nt. 

Hence t/gt is the projection of Ut onto span{5'(. i-g)}, which implies rank([/st) < s. 
With this in mind, if rank(t7gt) = then we have s^ < s. Otherwise, if rank(17st) < 
s^, we can conclude that s = rank([/gt) with the help of Lemma 13.11 and thereby 


6 
























s < s^. Based on these arguments, the following algorithm was proposed in [21] for 
finding a good upper bound for s. Meanwhile, a projection matrix onto the eigenspace 
span{S'(._]^.s)} is produced, which will play an important role in the resulting method 
given in next section. 

Algorithm 1. Input an increasing factor a > 1 and the size p of sample vectors. 
The function “Search” outputs si, an upper bound of the number of eigenvalues s 
inside F, and a projection matrix Ui £ onto span{S'(. i-g)}. 

Function [17i,si] = Search(A, i?, F, a,p) 

1. Pick Yp ~ N„xp(0,1) and compute U = ^ ^ {zB — A)~^BdzYp 

by the g-point Gauss-Legendre quadrature rule. 

2. Set So = |’-trace(F^*17)] and s* = min(max(p, sq), n). 

3. If s* > p 

4. Pick Y ~ N„x(s*-p)(0,1) and compute TJ = - ^-j= (f {zB — A)~^BdzY 

27rV —1 Jv 

by the g-point Gauss-Legendre quadrature rule. 

5. Augment G to C/ to form U = [U, U] £ 

6. Else 

7. Set s* = p. 

8. End 

9. Compute C/II = UxRx'. the rank-revealing QR decomposition |5] of U. 

10. Set Si = rank(i?). 

11. If Si < s*, stop. Otherwise, set p = si and s* = [asi]. Then go to Step 3. 

In practice, by (12. 7|) . we see that U (line 5) formed in the last iteration in Algo¬ 
rithm |T] is 



(3.2) 


where Yg* ^ N„xs* (0,1) with s* > si. According to Lemma ITTI we know that S "^Yg* 
is full-rank, then it follows from (13.211 that 


spanjUi} = span{[/} = span{S'(:^x)} 


(3.3) 


where X is an index set and its cardinality is si. It has been shown in previous section 
that > ^If^i, thus we can conclude that {1, 2,..., s} C X. 

4. Counting the eigenvalnes inside F. In this section, we present the result¬ 
ing method for counting the eigenvalues inside a given disk in the complex plane. 
Below is the second finding of our work. 

Theorem 4.1. Let Ui be the projection matrix computed by Algorithm]^ Let 


Compute U 2 by the q-point Gauss-Legendre quadrature rule, i.e., 




(4.1) 
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and define the Si x si matrix 


M = U 1 U 2 . (4.2) 

Then are the eigenvalues of M. 

Proof. Since Ui is computed by Algorithm [TJ by (13.31) . we can find a nonsingular 
matrix, say, W such that 


= UiW. 

Let T be the index set {1, 2,..., n} \I, then there exists a permutation matrix P such 
that 


Note that Z? is a diagonal matrix, thereby P*DP is also a diagonal matrix and can 
be written as 


where the diagonal entries of Di consist of S I. 

By dlH), (gl]) and g21), we have 

M = UfU2 = Uf{SDS-^)Ui 
= Uf{SPP*DPP*S-^)Ui 

= Uf{UiW)DiW-^ 

= WDiW-\ 


Since W is nonsingular, the matrices M and Di have the same eigenvalues, which are 
Note that {l,2,...,s} C I, therefore {Z?(i,i)}i=i are the eigenvalues of 

M. □ 

Theorem 14.11 tells us that can be obtained via computing the eigen¬ 

values of M (cf. (14.21) '). We have shown in Section 2 that the real parts of 
which correspond to the eigenvalues inside T, are larger than i, and the real parts of 

•1 ^ 

the rest diagonal entries of D are smaller than i. Motivated by these facts, we find 
the number s via counting the eigenvalues of M whose real parts are larger than i. 
Summarize this idea, below we give the complete algorithm for computing the number 
of eigenvalues inside T. 

Algorithm 2. Input an increasing factor a > 1 and the sizep of sample vectors. 
The function “Count_Eigs ” computes the number of eigenvalues of C3P that are 
located inside circle T. 

Function s = Count_Eigs(A, B, F, a,p) 

1. Call [17i,si] = Search(A, i?,r,a,p). 

2. Compute U 2 in (14.11) . and set M = 

3. Compute the eigenvalues of M, and set s to be the number of the computed 
eigenvalues whose real parts are larger than ^. 




As with other contour-integral based eigensolvers [B uni im na, the dominant 
work of Algorithm [5] is solving generalized shifted linear systems of form 


{zjB - A)Xj = BY, 


(4.3) 


which can be solved by any method of choice. Since the quadrature nodes Zj and 
the columns of the right-hand sides of (14.31) are independently, Algorithm [2] has good 
scalability in modern parallel architectures. 

The information about the number of eigenvalues inside the target region is crucial 
to the implementation of the contour-integral based eigensolvers. Algorithm [2] can 
integrate with this kind of eigensolvers [B [B HH [24] to provide them with this 
important information. For example, if we use the contour-integral based eigensolver 
proposed in [24] to compute the eigenvalues inside F, Ui and need to be computed 
in the first iteration and the second, respectively. As a result, to get s, the extra work 
is constructing the matrix M (cf. (14.21) 1 and then computing the eigenvalues of M. 
Once the number s is obtained, it can help to detect whether all eigenvalues inside 
F are found in the subsequent iterations and determine when to stop the iteration 
process. For clarity, we summarize the idea of integrating Algorithm [5] with the 
eigensolver proposed in [53] by the following algorithm. 

Algorithm 3. Input A,Bg an increasing factor a > 1, the size p 

of sample vectors, a circle F, a convergence tolerance e, and “maxAter” to control 
the maximum number of iterations. The function “Eigenpairs ” computes eigenpairs 
(Ai,Xi) of hl.l\) that satisfies 


Xi inside F and 


||Ax, - XiBx ,\\2 

||Axj|| 2 -b ||-Bxi||2 


(4.4) 


The results are stored in vector A and matrix X. 


Function [A, X\ = Eigenpairs(A, B, F, e, maxAter) 

1. Call [C/i,si] = SEARCH(A,i3,F, a,p). 

2. Compute U 2 in m, and let M = UfU 2 . Set s to be the number of the 
eigenvalues of M whose real parts are larger than i. 

3. For k = 2, ■ ■ ■ , maxAter 

4. Compute QR decompositions: Uk = UiRi and BUk = U 2 R 2 - 

5. Form A = UfAUi and B = U^BUi^ 

6. Solve the projected eigenproblem Ay = XBy of size si to obtain eigenpairs 

Set Xi = Uiyi,i = l,2,...,si. 

7. Set A = [ ] and AT = [ ]. 

8. For z = 1 : Si 

9. If (Ai,Xi) satisfies (14.4L then A = [A, Ai] and X = [X, x^j. 

10. End 

11. If there are s eigenpairs satisfying (14.4|) . stop. Otherwise, compute 

C/fc+i = ~ A)~^BdzUi by the g-point Gauss-Legendre 

quadrature rule: Uk+i ~ Uk+i = ^ ~ c)(-Zjl? — A)~^BUi. 

12. End 


Below we give some remarks on above algorithm. 

1. The first two steps can be viewed as using Algorithm [2] to determine the num¬ 
ber of eigenvalues inside F. In addition, a projection matrix is also produced, 
i.e., U 2 . 
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2. Since U 2 is a projection matrix onto span{S'(. it can be used to construct 
a projected eigenproblem to compute approximate eigenpairs. Therefore, in 
Step 3, the for-loop starts from k = 2. 

3. In Step 11, the number s computed in Step 2 helps to detect whether all s 
approximate eigenpairs satisfying the prescribed accuracy. If it is, we stop 
the iteration process. 

5. Numerical Experiments. In this section, we give some numerical experi¬ 
ments to illustrate the viability of our new method. All computations are carried out 
in Matlab version R2012b on a MacBook with an Intel Core i5 2.5 GHz processor 
and 8 GB RAM. 

In the experiments, as for computing the generalized shifted linear systems of the 
form (I4.3|l . we first use the Matlab function lu to compute the LU decomposition 
of ZjB — A, j = 1 , 2,..., g, and then perform the triangular substitutions to get the 
corresponding solutions. 

Experiment 5.1: Our method (Algorithmic) is motivated by two findings: (i) for 
matrix D (cf. (12.81) '). > 5 if correspond to the eigenvalues enclosed by 

r, and are less than ^ if Il(iy) correspond to the eigenvalues outside F; (ii) 

are the eigenvalues of M (cf. (14.2|) ') if they correspond to the eigenvalues within 
r. This experiment is devoted to illustrating these two findings. 

Let A = diag([0.1 : 0.1 : 0.8]), S = randn( 8 ), the matrices A and B are given by 


A = S'AS'-\ H = eye(8). 


Here diag, randn, and eye are Matlab commands. Obviously, for the problem 
under consideration, the eigenvalues are 0.1,0.2,..., 0.8. Let T be a circle with center 
at origin and radius p = 0.401. Suppose that we are interested in the number of 
eigenvalues inside T. Obviously, there are 4 eigenvalues are located inside T. Note 
that the eigenvalue 0.4 is located inside T and close to boundary of the disk surrounded 
by r. 

In this experiment, we select the number of quadrature points q = 32. According 
to ()2.8D . now D is given by 


32 



(5.1) 


where uji are the weights associated with quadrature nodes Zi- We take the size of 
sample vectors p = 6, thus the starting basis in function Search is Tg = randn(8, 6). 
Since the size of test problem is small and the number of columns of Ig is already larger 
than the number of eigenvalues inside F, we just run one iteration when preforming 
function SEARCH to get projection matrix Ui. As a result, the matrix M defined in 
g21) is of size 6x6. 

Since, in practical situations, we are only interested in the real parts of the di¬ 
agonal entries of D, in the second column in Table [5TT] we list that 

are computed by (15.11) . It can be seen that the corresponding to the 

eigenvalues inside F, are larger than 0.5. More precisely, = 1,2,3, corre¬ 

sponding to eigenvalues 0.1, 0.2, 0.3, respectively, approximate the theoretical value 1 
sufficiently. But 3fi[Zl(4q)] is about 0.8, this is because it corresponds to the eigenvalue 
0.4, which is close to the boundary of the target disk. On the other hand, the real 
parts of corresponding to the eigenvalues outside F, are less than 0.5 and 

very close to zero, as expected of our first finding. 
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Table 5.1 

The real parts of the diagonal entries of D and the ones of the eigenvalues of M. 




1 

2 

3 

4 

5 

6 

7 

8 


1.0000000000039 49 

1.0000000000000 00 

0.999999999999999 

0.8015817876596 01 

0.0000000025256 84 

0.0000000000043 79 

0.000000000000001 

-0.000000000000051 


SR[eig(A4)] 

1.0000000000039 65 

1.0000000000000 12 

0.999999999999999 

0.8015817876596 10 

0.0000000025256 20 

0.000000000004380 


Table 5.2 

Test problems from Matrix Market that are used in Experiment 5.2 


No. 

Matrix 

Size 

nnz 

Property 

cond 

1 

A: BFW398A 

B: BFW398B 

398 

398 

3678 

2910 

unsymmetric 
symmetric indefinite 

7.58 X lO'^ 
3.64 X lO^ 

2 

A: BFW782A 

B: BFW782B 

782 

782 

7514 

5982 

unsymmetric 
symmetric indefinite 

4.63 X lO"* 
3.05 X lOi 

3 

A: BCSSTK08 
B: BCSSTM08 

1074 

1074 

12960 

1074 

symmetric positive definite 
symmetric positive definite 

4.77 X 10' 
8.27 X 10® 

4 

A: BCSSTK27 
B: BCSSTM27 

1224 

1224 

28675 

28675 

symmetric positive definite 
symmetric indefinite 

7.71 X 10"^ 
1.14 X lO^o 

5 

A: PLAT1919 

B: PLSK1919 

1919 

1919 

17159 

4831 

symmetric indefinite 
skew symmetric 

1.40 X 10^® 
1.07 X lOi® 

6 

A: BCSSTK13 

B: BCSSTM13 

2003 

2003 

42943 

11973 

symmetric positive definite 
symmetric positive semi-definite 

4.57 X 10^® 

Inf 

7 

A: MHD4800A 
B: MHD4800B 

4800 

4800 

102252 

27520 

unsymmetric 
symmetric indefinite 

2.54 X 10®' 
1.03 X lOi"' 

8 

A: BCSSTK25 
B: BCSSTM25 

15439 

15439 

252241 

15439 

symmetric indefinite 
symmetric positive definite 

1.28 X 10^® 
6.06 X 10® 


The third column in Table IOI displays the real parts of the eigenvalues of M in 
descent order. The same digits of i)] and the ith largest 5R[eig(M)], i = 1,. .., 6 
are underlined. We can see that z = 1,..., 6, agree at least fourteen digits 

to their counterparts in the third column. Therefore the real parts of the eigenvalues 
of M are almost equivalent to q], * = 1,..., 6, which justifies the second finding 

of our work. 

Experiment 5.2: This experiment is devoted to testing the viability of our new 
method. The test matrices are available from the Matrix Market collection g]. They 
are the real-world problems from scientific and engineering applications. The descrip¬ 
tions of the matrices are presented in Table 15.21 where nnz denotes the number 
of non-zero entries and their condition numbers are computed by Matlab function 
condest. The test matrices of different problems vary in size, spectrum and property. 

In this experiment, we take the parameter q to be 16. The parameters c and p are 
the center and the radius of circle T, respectively. Note that the test problems 3, 4, 6, 
and 8 are Hermitian problems, which means their (finite) eigenvalues are real-valued. 
Due to this, we choose the circles with centers lying on the real line for these test 
problems. 

TABLE [5?3l Dresents the numerical comparisons, s is the actual number of eigenval¬ 
ues inside T. We first use the Matlab built-in function eig to compute all eigenvalues 
of the test problems, and then determine the values of s according to the coordinates 
of computed eigenvalues. We also present the estimation sq computed by the trace 
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Table 5.3 

Numerical comparison: s is the exact number of eigenvalues inside F, so is the estimate of s 
by using the trace formula and si is the upper bound computed by Algorithm^\ and Cont_Eigs is 
the result computed by our new method. 


No. 

(c, p) 

s 

so 

Si 

Cont.Eigs 

1 

((-6.0 X 10“) -1- V-l(2.0 X lO'^), 3.0 X lO'") 

120 

136 

177 

120 

2 

((-5.0 X 10®) -1- ^/^(1.0 X 10®), 2.0 X 10®) 

165 

161 

228 

165 

3 

(5.0 X 10®, 3.0 X 10®) 

178 

199 

190 

178 

4 

(5.0 X 10®, 3.0 X 10®) 

160 

134 

192 

160 

5 

(^/^(5.0 X 10-1), 1.5 X 10-1) 

301 

297 

397 

301 

6 

(6.0 X 10®, 3.5 X 10®) 

232 

229 

262 

232 

7 

((-5.0 X 10®) -1- ^/^(2.0 X 10®),4.0 X 10®) 

212 

545 

293 

212 

8 

(5.0 X 10®, 2.0 X 10®) 

1663 

1628 

1749 

1663 


formula (EU and the upper bound si computed by Algorithm [T] Cont_Eigs is the 
result computed by our new method. 

The results for all eight test problems are reported in Table I^TSI From these data, 
we see that the estimation sq computed by the trace formula (ED always provide good 
estimation to s for all test problems except for test problem 7, due to its ill-condition. 
Si always gives a good upper bound for s. It is remarkable that the result computed 
by our new method is the same with the exact number s for each test problem, even 
though the test problem 7. Therefore, our new method is numerically efficient and 
reliable. 


6. Conclusion. In this work, we develop an approach for counting the eigenval¬ 
ues of mD inside a given disk in the complex plane. The new method is a contour- 
integral based method and motivated by two findings. The computational advantage 
of the new method is that it is easily parallelizable. Its another promising feature is 
that it can integrate with the recently proposed contour-integral based eigensolvers to 
provide them the information of the number of eigenvalues inside the target region. 
How to adapt the resulting method to the nonlinear problems will be our future work. 
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